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Abstract 

r~| ■ We study the thermodynamic Casimir force for films in the three-dimensional Ising 

^ , universality class with symmetry breaking boundary conditions. To this end we simulate 

Ch ! the improved Blume-Capel model on the simple cubic lattice. We study the two cases 

++, where all spins at the boundary are fixed to +1 and +— , where the spins at one 
C^ ■ boundary are fixed to +1 while those at the other boundary are fixed to —1. An important 

^ ■ issue in analyzing Monte Carlo and experimental data are corrections to scaling. Since 

j^ , we simulate an improved model, leading corrections to scaling, which are proportional 

H ! to Lq'^ , where Lq is the thickness of the film and to ~ 0.8, can be ignored. This allows 

_L I us to focus on corrections to scaling that are caused by the boundary conditions. The 

^ ■ analysis of our data shows that these corrections can be accounted for by an effective 

^ ■ thickness -^o,e// = -^o + Ls- Studying the correlation length of the films, the energy per 

area, the magnetization profile and the thermodynamic Casimir force at the bulk critical 
^vq ! point we find Lg = 1-9(1) for our model and the boundary conditions discussed here. 

^ I Using this result for Lg we find a nice collapse of the finite size scaling curves obtained 

for the thicknesses Lq = 8.5, 16.5 and 32.5 for the full range of temperatures that we 






C^ ■ consider. We compare our results for the finite size scaling functions 0_|__|_ and 9^ of the 



^ 



X 



thermodynamic Casimir force with those obtained in a previous Monte Carlo study, by 
the de Gennes-Fisher local-functional method, field theoretic methods and an experiment 



o 

C^ [ with a classical binary liquid mixture. 
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I. INTRODUCTION 

In the thermodynamic hmit, in the neighborhood of a second order phase transi- 
tion the correlation length ^ that is the characteristic length of thermal fluctuations 
diverges following a power law 

^ = ^o,±\trx{i+b±\tf+ct+...) , (1) 

where t = (T ~ Tc)/Tc is the reduced temperature and ^o,± is the amplitude of the 
correlation length in the low (— ) and the high (+) temperature phase, respectively. 
Using this notation, we assume that the high temperature phase is characterized 
by disorder and the low temperature one by order. The power law ([T]) is subject 
to confluent corrections, such as &±|t|^, and non-confluent ones such as ct. Critical 
exponents like u and ratios of amplitudes such as ^o,+/^o,- are universal. This 
means that they assume exactly the same value for any system within a given 
universality class. Also correction exponents like 6 = uu and ratios of correction 
amplitudes as b^/b^ are universal. For the three-dimensional Ising universality, 
which is considered here and other three-dimensional universality classes like the 
XY or the Heisenberg universality class, 6 ^ 0.5. For reviews on critical phenomena 
and the Renormalization Group (RG) see e.g. \Mi^- 

In 1978 Fisher and de Gennes J5| realized that when thermal fluctuations are 
restricted by a container, a force acts on its walls. Since this effect is analogous to 
the Casimir effect, where the restriction of quantum fluctuations induces a force, it 
is called "thermodynamic" Casimir effect. Since thermal fluctuations only extend 
to large scales in the neighborhood of continuous phase transitions it is also called 
"critical" Casimir effect. Recently this force could be detected for various experi- 
mental systems and quantitative predictions could be obtained from Monte Carlo 
simulations of spin models |6|. 

Here we study the thermodynamic Casimir force for the film geometry. From a 
thermodynamic point of view, the thermodynamic Casimir force per area is given 
by 
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where Lq is the thickness of the film and fex = ffUm — Lofbuik is the excess free 
energy per area of the film, where ffUm is the free energy per area of the film and 
fbuik the free energy density of the bulk system. The thermodynamic Casimir force 
per area follows the finite size scaling law 
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ksTL,' 9it[Lo/U+Y^n , (3) 



see e.g. ref. ^]. The finite size scaling function 6{x) depends on the universality 
class of the bulk phase transition, the geometry of the finite system and the surface 
universality classes of the boundary conditions that are applied. For reviews of 
surface critical phenomena see |8l-ll0|. Similar to the power law ([1]), finite size 



scaling equations such as eq. (E]) are subject to corrections to scaling. In the generic 
case one expects that leading corrections are oc Lq'^ (Ref. [ll|), where u = 0.832(6) 
(Ref. 12|) for the three-dimensional Ising universality class. Furthermore one 



expects corrections that are caused by the boundaries. We shall give a more detailed 
discussion of corrections to scaling below in section IIVI 

Here we compute finite size scaling functions 9 of the thermodynamic Casimir 
force for the three-dimensional Ising universality class and symmetry breaking 
boundary conditions. Experimentally this situation is realized for example by a film 
of a classical binary liquid mixture. Typically, the surface is more attractive for one 
of the two components of the mixture, breaking the symmetry at the boundary. In 
the Ising model this can be described by an external field that acts on the spins at 
the surface of the lattice. Following the classification of surface critical phenomena 
such surfaces belong to the normal surface universality class, which is equivalent to 
the extraordinary surface universality class [l3|. In recent experiments on colloidal 
particles immersed in a binary mixture of fiuids [l4l], the authors have demonstrated 
that the adsorption strength can be varied continuously by a chemical modification 
of the surfaces. In particular the situation of effectively equal adsorption strengths 
for the two fiuids can be reached. For sufficiently small ordering interaction at the 
surface, this corresponds to the ordinary surface universality class. Hence these 
experiments open the way to study the crossover between different surface univer- 
sality classes. For a recent theoretical discussion of the crossover behaviors of the 
thermodynamic Casimir force see [l5| and references therein. Here we shall not 
study such crossover behaviors and restrict ourself to compute finite size scaling 
functions for the normal or extraordinary universality class. Note that the breaking 
of the effective symmetry between the components of the fluid, or the breaking of 
the Z2 symmetry between + and — spins at the surface in the Ising model, con- 
stitutes a relevant perturbation at the ordinary flxed point |8l-ll0|. Therefore, even 
for a small breaking of the symmetry, for sufficiently large distances, which means 
in our context a large thickness of the fllm, the physics in the neighborhood of the 
critical point is governed by the normal or extraordinary universality class. 

Since a fllm has two surfaces, we can distinguish the two principal cases: Firstly 
both boundaries attract positive spins, denoted by ++ in the following, and sec- 
ondly one boundary attracts positive spins, while the other attracts negative spins, 

denoted by H — in the following. Note that by symmetry and — \- boundary 

conditions are equivalent to ++ and H — boundary conditions, respectively. 

In previous Monte Carlo studies la, ll7| the spin-1/2 Ising model has been sim- 
ulated. Computing flnite size scaling functions from numerical data obtained for 
flnite thicknesses Lq, corrections to scaling are a major obstacle. The results for 
^++ and 6'+_ given by [ig, [iTj depend quite strongly on the ansatz that is chosen 



for the corrections. Here we shall study the improved Blume-Capel model on the 
simple cubic lattice. The Blume-Capel model is a generalization of the Ising model. 
In addition to ±1, as in the Ising model, the spin might assume the value 0. The 
parameter D of the model controls the relative weight of and ±1. For a precise 
deflnition see section [TTl below. Improved means that the amplitude of corrections 
oc Lq^ vanishes or in practice is very small compared with the spin-1/2 Ising model. 
Studying thin fllms this is a quite useful property, since the boundary conditions 
cause corrections that are oc L^^ as we shall discuss below. Fitting numerical data, 
it is quite difficult to disentangle corrections that have similar exponents. Avoiding 
this problem we are able to compute the finite size scaling functions ^+_|_ and ^+_ 
with a small and, as we shall argue, reliable error estimate. Reliable numerical cal- 



culations are important, since field theoretic methods do not provide quantitatively 
accurate results for the scaling functions 6^++ and 6'+_ as we shall see below. Re- 
cently the scaling function ^_|_+ has been computed by using the de Gennes-Fisher 



local-functional method [19[. We find a rather good agreement with our result. 

The outline of the paper is the following. First we define the model and the 
observables that we have studied. Then we discuss finite size scaling and corrections 
to finite size scaling. Next we exploit the relation of the spectrum of the transfer 
matrix and the thermodynamic Casimir force. Then we discuss the Monte Carlo 
algorithms that we have used. We analyze our data obtained from simulations at 
the critical point of the bulk system. This way we obtain accurate results for the 
Casimir amplitudes and for Lg that characterizes the corrections to scaling caused by 
the boundary conditions. Next we have simulated in a large range of temperatures 
around the bulk critical point. Based on these simulations we obtain the finite size 
scaling functions 6*++ and ^+_ of the thermodynamic Casimir force. In addition 
we compute the finite size scaling functions of the correlation length of the films. 
Finally we compare our results with those obtained by field theoretic methods, the 
local-functional method, previous Monte Carlo studies of the Ising model and an 
experiment on a classical binary liquid mixture. 

II. MODEL 

We study the Blume-Capel model on the simple cubic lattice. It is defined by 
the reduced Hamiltonian 

<xy> X 

where the spin might assume the values s^ G { — 1,0,1}. x = (xo,a;i,X2) denotes 
a site on the simple cubic lattice, where Xj G {l,2,...,Lj} and < xy > denotes a 
pair of nearest neighbors on the lattice. The inverse temperature is denoted by 
P = l/ksT. The partition function is given hj Z = Ylis} ^^pI""-^)) where the sum 
runs over all spin configurations. The parameter D controls the density of vacancies 
Sx = 0. In the limit D — )■ — cxd vacancies are completely suppressed and hence the 
spin- 1/2 Ising model is recovered. 

In d > 1 dimensions the model undergoes a continuous phase transition for 
— oo < D < Dfri at a Pc that depends on D. For D > Dtri the model undergoes 
a first order phase transition. The authors of [20] give for the three-dimensional 
simple cubic lattice D^ri = 2.0313(4). 

Numerically, using Monte Carlo simulations it has been shown that there is a 
point {D*, Pc{D*)) on the line of second order phase transitions, where the amplitude 
of leading corrections to scaling vanishes. Our recent estimate is D* = 0.656(20) 

I — H I — —I ^ ^ ^ 

(Ref. [12[). In [12[ we have simulated the model at D = 0.655 close to (3c on 
lattices of a linear size up to L = 360. From a standard finite size scaling anal- 
ysis of phenomenological couplings like the Binder cumulant we find /3c(0.655) = 
0.387721735(25). Furthermore the amplitude of leading corrections to scaling is at 
least by a factor of 30 smaller than for the spin-1/2 Ising model. 

In [2l| we have simulated the Blume-Capel model at D = 0.655 in the high 
temperature phase on lattices of the size L^ with periodic boundary conditions in 



all directions and L ^ 10^ for 201 values of /9. We have measured the second moment 
correlation length ^2nd that we shall define below. The simulation at /3 = 0.3872, 
which was our closest to /3c, yielded ^2nd = 26.698(7). Fitting these data for ^2nd 
with ansatze obtained by truncating the sequence of correction terms at various 
order we arrive at 

6nd,o,+ = 0.2282(2) - 1.8 x (z/ - 0.63002) + 250 x {(3c - 0.387721735) 

using t = f3c — (3 as definition of the reduced temperature. (5) 



In these fits we have fixed u = 0.63002 and /Sc = 0.387721735 (Ref. [jj). We have 
redone the fits with slightly shifted values of z/ and /3c to determine the dependence 
of ^2ndfl,+ on these input parameters. For simplicity we shall use t = /9c — /9 as 
reduced temperature also in the following. 

In the high temperature phase there is little difference between ^2nd and the 
exponential correlation length C,exp which is defined by the asymptotic decay of the 



two-point correlation function. Following [22 



lim|^ = 1.000200(3) (6) 

for the thermodynamic limit of the three-dimensional system. This means that at 
the level of our accuracy we can ignore this difference. Note that in the following 
^0 always refers to ^2nd,o,+ , eq. ©. 



A. Film geometry and boundary conditions 

In the present work we study the thermodynamic Casimir effect for systems with 
film geometry. In the ideal case this means that the system has a finite thickness Lq, 
while in the other two directions the thermodynamic limit Li, L2 — !■ oo is taken. In 
our Monte Carlo simulations we shall study lattices with Lq <^ Li, L2 and periodic 
boundary conditions in the 1 and 2 directions. Throughout we shall simulate lattices 
with Li = L2 = L. 

In the direction we take symmetry breaking boundary conditions. In the 
reduced Hamiltonian of the Blume-Capel model these can be implemented by 

H = -(3 ^ Sa:Sy + D^sl - hi ^ Sx - h2 ^ s^ , (7) 

<xy> X 2:0=0,2:1, X2 xo=Lo+l,xi,X2 

where hi,h2 7^ break the symmetry at the surfaces that we have put on xq = 
and Xq = Lq + 1. Hence Lq gives the number of layers in the interior of the film. 

In our Monte Carlo simulations we consider the limit of infinitely strong surface 
fields hi and /12, which means that the spins at the surface are fixed to either —1 
or 1, depending on the signs of hi and /i2- Therefore we have implemented in our 
simulation code ++ boundary conditions by setting s^^ = 1 for all x with xq = or 
Xq = Lq + 1 and H — boundary conditions by setting Sx = I for all x with xq = 
and Sx = —^ for all x with xq = Lq + 1. Alternatively, these fixed spins could be 
interpreted as finite surface fields with \hi\ = I/12I = (3 acting on the spins at Xq = 1 
and Xq = Lq, respectively. 



III. OBSERVABLES 

A. Internal energy and free energy 

The reduced free energy per area is defined by 

/ = --i-lnZ. (8) 

Tliis means that compared with the free energy per area /, a factor ksT is skipped. 
Correspondingly we define the energy per area as the derivative of minus the 
reduced free energy per area with respect to /3: 

1 dlnZ 1 / v^ \ 

It is straight forward to determine E in Monte Carlo simulations. From the defini- 
tion of E follows 

/(/3) = /(/3o) - / dpE0) . (10) 

B. The magnetization profile of films 

The film is invariant under translations in the 1 and 2 direction of the lattice. 
Therefore the magnetization only depends on xq and we can average over xi and 
X2: 

"^(^0) = j2 ^{sx) ■ (11) 

Since the film is symmetric for ++ boundary conditions and anti-symmetric for H — 
boundary conditions under reflections at the middle of the film, itlIxq) = 'm{Lo — 
Xo + l) for -|--|- boundary conditions and m(a;o) = —mlLo — XQ + l) for H — boundary 
conditions. 



C. Second moment correlation length of the films 

We have measured the second moment correlation length of the films in the 1 and 
2 direction of the lattice. To this end we have computed the connected correlation 
function of the Fourier transformed field 

Gih,k2) = {mh,k2)\^) - 5^k,M),mLoL^^^ (12) 

where m is the magnetization and the Fourier transformed field 

ip{ki,k2) = ^===^expii ^^ ^-^jsx- (13) 



For large L and small fci, ^2? the correlation function behaves as 

C 

^^^"^'^ ^ 4sin2(7rfci/L) + ^sin\nh/L) + ^2"', ' ^^^^ 

The second moment correlation length ^2nd can now be evaluated by computing 
G(fci, ^2) for two values of (/ci, ^2) and solving eq. flT4|) with respect to ^|„^. In the 
limit L —7- 00 all choices of (/ci,/c2) lead to the same result for Q^^. However, for 
finite L the deviations from this limit increase with increasing values of ki and ^2- 
Therefore, for H — boundary conditions, we have computed the correlation function 
at (/ci,/c2) = (0,0) and (1,0). One gets 



^2 _ G(0,0)/G(1,0) -1 
4sin2(7r/L) 






In the simulation we have also measured (5(0,1) and have averaged G(1,0) and 
G(0, 1) to reduce the statistical error. 

In contrast to H — boundary conditions, for ++ boundary conditions there is 
a finite magnetization at any finite temperature. In order to avoid the technical 
complication of subtracting the magnetization squared required for (fci, k'l) = (0, 0), 
eq. ( 1T^ . we have used {ki, ^2) = (1, 0) and (1, 1) to determine the second moment 
correlation length 

.2 G(1,0)-G(1,1) 

'"'^ [2^(1,1) -G(1,0)] 4 sin2(7r/L) ' ^ ' 

In the simulations below we have chosen the lattice size L such that the limit L — > 00 
is well approximated. Hence ^2nd is a function of the parameters (3 and D of the 
model and the thickness Lq of the film. 



IV. FINITE SIZE SCALING 

The reduced excess free energy of the film behaves as 

fe.{Lo,t) = ffam{Lo,t) - Lohuikit) ^ L,''+^h{t[L,/ioY'') , (17) 

where ffiim{LQ,t) is the reduced free energy per area of the film, fhuikif) the re- 
duced free energy density of the bulk system, /i(t[Lo/^o]^'''^) is the universal finite 
size scaling function of the excess free energy and rf = 3 is the dimension of the 
bulk system. Here and in the following ^0 is the amplitude of the second moment 
correlation length of the bulk system in the high temperature phase. 

Inserting the finite size scaling ansatz (fT7|) for the excess free energy into (|2]) one 
gets 

^ Casimir — n^is-' 
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-{d - l)hit[Lo/^o]'^n + lt[Lo/^oY^^h\t[Lo/^oY^n 



where 

e{x) = {d-l)h{x)--h'{x) (19) 

V 

is the finite size scahng function of the thermodynamic Casimir force and x = 
t[Lo/^o]^'^''- This relation is well known and can be found e.g. in [7|. 

Following the discussion in section III B of ref . [ll| , taking into account leading 
corrections to scaling one gets 

fe,{Lo,t) = 1^"+^ h{x,a{D)L^'') = L^'^^Kx) x (1 + a{D)c{x)L^^ + ...) (20) 

and correspondingly for the thermodynamic Casimir force per area 

Fcas^m^r = kBTL^" ^(x, a{D)L^^) = kBTL^'' Oix) X {l+a{D)d{x)Lo'^ + ...) , (21) 

where we have performed the Taylor expansion of h and 6 in their second argument 



to leading order. The authors of [16|, |17[ arrive at a similar expression as eq. fl2T|) . 



Fitting their data, obtained for the Ising model, they have approximated the func- 
tion d{x) by a constant. For the improved model that we study here a{D) ^ 
holds, which simplifies the analysis of our data. 

The exponent of the leading correction to scaling takes the value u = 0.832(6) 
(Ref. [12] )■ Furthermore there are subleading corrections. Among these, the leading 
ones come with the exponents u' = 1.67(11) (Ref. |23|) and due to the breaking 
of rotational symmetry by the lattice to" ~ 2 (Ref. [2J]). At the level of accuracy 
of our data, we can not resolve the individual subleading corrections. In order to 
get some estimate of the effect of these corrections on our final results, we have 
included a term cLq"^ into the ansatze (I37f39f45ll48l) below. 

A discussion of corrections caused by the boundaries is given in section V A of 
ref. [ll|. Corrections might arise from irrelevant surface scaling fields. Furthermore 
Capehart and Fisher [25| have argued that there is an arbitrariness in the definition 
of the thickness of the film leading to corrections oc Lq^. These two arguments 
might be actually unified: In a real-space Renormalization Group treatment of 
surface critical phenomena one splits the reduced Hamiltonian into a bulk and a 
surface part. In the neighborhood of the critical point, one might expand the bulk 
and the surface part of the reduced Hamiltonian into so called scaling fields. The 
basic idea is that splitting the reduced Hamiltonian into a bulk and a surface part 
is a priori quite ad hoc. Roughly speaking, one might put the contribution for 
(1 — ls)/2 < Xo < i^o + (1 + ^s)/2 of eq. (jl]) into the bulk part and the remainder 
into the surface part. This way, the amplitudes of the surface scaling fields become 
functions of Ig- Here we do not elaborate what sense can be given to non-integer 
values of Ig- The amplitude of the leading irrelevant surface scaling field, viewed 
as a function of Ig, might have a zero that we shall call Lg in the following. Then 
this surface scaling field has the RG exponent yg = —Ug = — 1. If there is only one 
surface scaling field with the RG exponent yg = —1, corrections oc Lq^ can hence 
be eliminated by replacing Lq by Lo^eff = Lq + Lg in finite size scaling laws. 

For the ordinary surface universality class, the problem of corrections has been 



worked out in some detail. A field theoretical calculation [26[ predicts a single ir- 
relevant scaling field with the RG exponent yg = —1. These corrections to scaling 
are related with the extrapolation length, which was introduced in the context of 



mean-field theory; See the review [8[. It is given by the zero of the extrapolated 



magnetization profile. The authors of [27] have employed the concept of the ex- 
trapolation length in their Monte Carlo study of the magnetization profile of the 
three-dimensional Ising model on the simple cubic lattice with free boundary condi- 
tions, which belong to the ordinary surface universality class. They have simulated 
various values of the ratio w of the surface and the bulk coupling. They find that 
the data for different values of w only fall nicely on a single scaling curve, when the 
extrapolation length that depends on w is properly taken into account. Finally we 
like to mention that there had been attempts to eliminate corrections due to the 
surface by a proper choice of w |28| . 

It is beyond the scope of the present manuscript to check whether the result 
of the field theoretical calculation [26] carries over to the extraordinary surface 
universality class, which is relevant for the present study. Our working hypothesis 
is that there is only a single irrelevant surface scaling field with the RG exponent 
Us = —1 which can be accounted for by an effective thickness L^^eff of the film. 
Furthermore we assume that there are no other irrelevant surface scaling fields with 
ys ~ —2. The analysis of our precise numerical data for various quantities provides 
a quite non-trivial challenge of this hypothesis. 

Finally let us spell out how the effective thickness -Z^o.e// enters into finite size 
scaling laws. For the thermodynamic Casimir force one gets 

Fcasimir = ^bT L^^^^ ^0{t[LQ^ef f / ^o] ) (22) 

where both the prefactor L^ as well as the scaling variable x = t[Lo/^o]^'''^ are 
replaced by L^eff and x = t[Lo,e///^o]^^'^) respectively. 

We also study the finite size scaling behavior of the second moment correlation 
length of the film. Taking into account boundary corrections we get 

^2nd,film = Lo^effX{t[Lo^eff/^o] ) • (23) 

The magnetization profile at the bulk critical point behaves as 

m{xo) = c L'^Jj) ipiz/Lo^eff) , (24) 

where z = xo — Lq/2 — 1/2 gives the distance from the middle of the film and c is 
a model specific constant that could be fixed by the behavior of the magnetization 
or the magnetic susceptibility in the thermodynamic limit. From scaling relations 
it follows that /3/z/ = (1 -|- r])/2, where rj = 0.03627(10) for the three-dimensional 
Ising universality class [l2[. Note that the scaling function 'ip{z/Lo^eff) diverges as 
z/Lo^eff ~^ il/2! since the magnetization in the neighborhood of the boundary 
stays finite as Lq — )■ oo for the boundary conditions studied here. 



A. Thermodynamic Casimir force and the transfermatrix 

The partition function of the system with fixed boundary conditions can be ex- 
pressed in terms of the eigenvalues of the transfermatrix and the overlap of the 
eigenvectors with the boundary states. Let us consider a lattice of the size LqX L"^, 
where L is large compared with the bulk correlation length but still finite. We con- 
sider the transfermatrix T that acts on vectors that are build on the configurations 



living on L^ slices. We denote the eigenvalues of T by Aq and the corresponding 
eigenvector by |a), where a = 0, 1, 2, ..., a^ax- The eigenvalues are ordered such 
that Aq, > A/3 for a < /3. Note that T commutes with translations, rotations, re- 
flections and with the change of the sign of all spins in a slice. Therefore the states 
I a) can be classified according to their momentum, the angular momentum, their 
parity and their behavior under sign-change of the spins. Note that on the lattice, 
only a sub-group of the symmetries of the continuum is realized. For a detailed 
discussion of the implications of this fact see for example section 3.2 of 29j, where 
the spectrum of the Ising gauge model in 2 + 1 dimensions had been studied. 
Now we can write the partition function of the system with fixed boundaries as 

Zt,,b, = J2^'o.{bi\a){h\a) , (25) 

a 

where I = Lq + 1 for our definition of the thickness Lq. The boundary states 61,2 
can be either + or — here. Note that these boundary states are invariant under 
all symmetries discussed above except for the sign-change of the spins. Therefore 
only states |a) with zero momentum, zero angular momentum and even parity have 
a non- vanishing overlap {b\a). Now we can compute the thermodynamic Casimir 
force per area starting from eq. fl2Sl) 

1 Id 
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(26) 



L^ EJAJAo)' {h\a){h2\a) 

where Aq is the largest eigenvalue. Introducing the inverse correlation lengths 
l/ia = ma = - ln(A„/Ao) we get 

1 p ^ ^ ^^m„exp(-m^/) {bi\a){b2\a) ,^7) 

ksT "'''"'''' L2 ^^exp(-m„/) {bi\a){b2\a) 

This equation proves that for 61 = 62 the thermodynamic Casimir force takes neg- 
ative values. In the high temperature phase, in the zero momentum sector, the 
second largest eigenvalue Ai is well separated from larger eigenvalues. Therefore 
the behavior of the thermodynamic Casimir force for / ^ ^1 = ^ = 1/m, which 
corresponds to large values of the scaling variable x, is given by 



9{ml) 
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ksT "'^^^""'^ (6i|0)(62|0)+exp(-m/)(6i|l)(62|l) 



.3,3 _. ^n 1 (&i|l)(&2|l) 



.-rfexp(-^0;;^|^^^^. (28) 

The finite size scaling behavior f lTSjl of the thermodynamic Casimir force implies 

that 

Ci'^ = 4M (29) 

^ ' mL (6|0) ^ ' 

has a finite scaling limit. The state |0) is symmetric under the global transformation 
Sx — )■ —Sx for all X in a slice. Instead, |1) is anti-symmetric and therefore C = 
C{+) = -C(-). It follows 

e++{ml) = -e+_{ml) = -C^ mH^exp{-ml) (30) 



for sufficiently large values of ml. Since x = t[l/^QY^'^ ~ [mlY^^ it follows 

e^^{x) = -e+-(x) = -C^x^'^expi-x") (31) 

for sufficiently large values of x. In the low temperature phase, the situation is more 
complicated. Also here, for finite L the state |0) is symmetric under Sx -^ —Sx, while 
|1) is anti-symmetric. The corresponding correlation length ^j = — l/ln(Ai/Ao) is 
the so called tunneling correlation length. It diverges as ^t oc exp(crL^) in the 
limit L — !• oo, where a is the interface tension. It is characteristic for the low 
temperature phase, and a consequence of spontaneous symmetry breaking that pairs 
of eigenvalues, where one is symmetric and the other anti-symmetric under Sx — )■ 
—Sx, become degenerate in the limit L — t- oo. The bulk correlation length in the low 
temperature phase is given by ^ = — limi-^oD 1/ ln(A2/Ao) = — limi^oo 1/ ln(A3/Ao). 
Taking into account the states a = 0, 1, 2 and 3 we get 

1 ^ 1 m2exp(-m20(&i|2)(62|2)+m3exp(-m3/)(6i|3)(62|3) 
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ksT '^'^^^^^ L2 (6i|0)(62|0)+exp(-m,/)(6i|l)(62|l) 

(32) 
where we have skipped the contribution of a = 1 in the numerator, since rrit vanishes 
in the limit L — )■ oo. Furthermore, we have skipped the contributions of a = 2 and 
3 in the denominator, since for 777-2/, f^s^ ^ 1 they are small compared with those 
of a = and 1. For H — boundary conditions (-|-|a)(— |q;) is positive for states 
that are symmetric and negative for states that are anti-symmetric under the spin- 
flip. Therefore both in the numerator and the denominator there is a cancellation 
between the two terms. Extracting useful information from eq. fl32|) would require 
detail knowledge of the approach of rrit, 7772, 7^73 and the overlap amplitudes to the 
limit L — )■ 00. 

On the other hand for ++ boundary conditions (+|Q;)(+|a) is positive for any 
a. Therefore in eq. ( l32l) the two terms in the numerator and the denominator add 
up. In the limit L — )■ 00, where rrit = and 7r7 = 7r72 = 7773 we get a result analogous 
to eq. (130!) . We only have to notice that in the definition of the scaling variable x 
the amplitude ^0,+ of the correlation length in the high temperature phase enters. 
Therefore taking into account the universal amplitude ratio ^o,+/^o,- = 1.901(14) 
(Ref. [2l|) for the exponential correlation length we get 

e++{x) = -C^[-1.901{U)xf'exp{-[-1.901{U)xY) (33) 

for sufficiently small values of x in the low temperature phase. For a discussion of 
the spectrum and the symmetry properties of the eigenvectors of the transfermatrix 
see e.g. |30[. Eqs. fl3H [33]) had been derived before by using the de Gennes-Fisher 



local- functional method, see eq. (6) of ref. [19| . Exact results for the Ising strip [31 
and mean-field theory |32[ confirm the exponential decay of 6^+{x) for large |a;|. 



V. MONTE CARLO ALGORITHMS 

A. ++ boundary conditions 

In the case of ++ boundary conditions we have used a hybrid of a cluster update 



and a local heat bath algorithm [33|. The cluster algorithm can only change the 



sign of the spins. Therefore local heat bath updates are needed to get an ergodic 
algorithm. For the cluster algorithm, we have used the same probability to freeze 
or delete a link < xy > as it is used in the original Swendsen-Wang [3^ algorithm: 

Pd{s:,:Sy) =mm[l,exp{-2f3s:rSy)] . (34) 

Links are deleted with the probability Pd{sxSy), otherwise they are frozen. A cluster 
is a set of sites that is connected by frozen links. In the following we mean by 
"flipping a cluster" that the sign of all spins Sx, where the site x belongs to the 
cluster, is changed ("flipped"). In one step of the Swendsen-Wang cluster algorithm, 
the lattice is completely decomposed into clusters. A cluster is then flipped with 
the probability 1/2. In contrast, in the case of the Wolff single cluster algorithm 
35|, one site of the lattice is chosen randomly. Then only the cluster that contains 



this site is constructed. This cluster is flipped with probability 1. Here we have 
to deal with the boundaries. For links < xy >, where either x or y belongs to 
the boundary we shall apply the same freeze or delete probability fjM|) as for links 
< xy >, where none of the two sites belongs to the boundary. Since spins on the 
boundary are fixed to one, clusters that contain sites on the boundary can not be 
flipped. Motivated by this fact, we have flipped all clusters with probability one 
that do not include sites on the boundary. In practice this is done in the following 
way: First we compute all clusters that include sites on the boundary. Then all 
spins on sites that do not belong to these clusters are flipped. 

With the local heat bath algorithm we run through the lattice in typewriter 
fashion. Running through the lattice once is called one "sweep" in the following. 
One cycle of the hybrid algorithm is composed of two sweeps of the local heat bath 
algorithm followed by one cluster update as discussed above. At the bulk critical 
point the integrated autocorrelation time of the energy is Tint,E ~ 3 in units of 
update cycles for a lattice of the size_Lo = 32, Li = L2 = 128. The integrated 
autocorrelation times for G{1, 0) and G{1, 1) are smaller. 



B. H — boundary conditions 

We could not use the program written for the ++ boundary conditions for the 
H — boundary conditions, since it relies on the fact that all spins that belong to 
clusters that include sites on the boundary are equal to +1. For simplicity we 
therefore have used a local Me trop olis algorithm that was implemented by using 
the multispin coding technique |36| . Details of our implementation can be found in 
[l2|. In [l2| we have found a performance gain of our Metropolis update using the 
multispin coding technique of about a factor of ten compared with the heat bath 
algorithm, implemented in a standard way. 

Likely, for small values of Lq the local Metropolis algorithm implemented by 
using the multispin coding technique outperforms the hybrid of local heat bath and 
cluster algorithm in the case of ++ boundary conditions. For lack of time we did 
not check this. 

In the low temperature phase, for H — boundary conditions rather large auto- 
correlations arise. These are due to fluctuations of the interface between the + and 
the — phase. As discussed in |37| standard cluster algorithms are not suitable to 



overcome this problem. Unfortunately, the algorithm discussed in 37| only works 
well in the Ising limit. 

In all our simulations we have used the SIMD-oriented Fast Mersenne Twister 



algorithm [38| as random number generator. 



VI. SIMULATIONS AT THE BULK CRITICAL POINT 

Here we focus on the finite size scaling behavior of various quantities at the 
bulk critical point. This way we accurately compute L^, which characterizes the 
corrections caused by the boundary conditions. To this end we have performed two 
sets of simulations. First we have simulated films of the size Lq x L^ to determine 
the second moment correlation length in 1 and 2 directions, the energy per area of 
the films and the magnetization profile. Then we computed the differences 

A/(Lo, /?,) = /(Lo + 1/2, /?,) - /(Lo - 1/2, /3,) (35) 

of free energies per area, where Lq + 1/2 and Lq — 1/2 assume integer values. To this 
end, we have simulated a lattice with Lq — 1/2 complete layers and one incomplete 
layer. A/(Lo, /3c) is then given by the free energy required to add a single site to 
this incomplete layer. For details of the method see |39| . 



A. Correlation length and energy per area at the bulk critical point 

For both H — and ++ boundary conditions we have simulated lattices of the 
thicknesses Lq = 6, 7, 8, ..., 26, 28, 30, 32. Throughout we have used L = ALq. At 
the bulk critical point, the correlation length of films with ++ boundary conditions 
is C,2nd ~ 0.13Lo and for H — boundary conditions ^2nd ~ 0.21Lo, as we shall see 
below. Therefore this choice of L is sufficient to get a good approximation of the 
limit L — )■ cxD. Throughout we have performed 100 000 000 update cycles for ++ 
boundary conditions and 64 x 5 000 000 measurements for H — boundary conditions. 
In the case of H — boundary conditions up to 18 Metropolis sweeps were performed 
for each measurement. In total the simulations took one year and 1.5 years on one 
core of a Quad-Core AMD Opteron(tm) Processor 2378 running at 2.4 GHz for ++ 
and H — boundary conditions, respectively. 

We have ffited the second moment correlation length at the critical point of the 
bulk system with the ansatz 

^2nd = c (Lo + Ls) (36) 

and to check for the possible effect of subleading corrections 

^2nd = c (Lo + L,) X (1 + 6 (Lo + L,y^) . (37) 

For ++ boundary conditions, ffiting with ansatz (!36|) we get for Lq^iiu = 12 the 
resuhs c = 0.1303(2), L, = 1.89(3) and xVd-o.f.= 0.83. In this ffi we have taken 
all data with Lq > Lo,mm into account. Using instead the ansatz (I57|) we get for 
Lo,min = 6 the results c = 0.1303(2), L, = 1.89(4) and xVd.o.f.= 0.94. 



For H — boundary conditions, fitting with ansatz fl5^ we get for Lo,mm = 14 
tlie results c = 0.2111(3), L, = 2.01(3) and xVd.o.f.= 1.97. Using instead tlie 
ansatz (l37j) we get for L^rnin = 8 the resuhs c = 0.2119(4), Lg = 1.81(6) and 
X^/d.o.f.= 2.08. In both cases, the x^/d.o.f. does not further decrease with increas- 
ing Lo,mm- 

We conclude that the results obtained for Lg for the ++ and the H — boundary- 
conditions are both consistent with Lg ^ 1.9. We have checked that the error of (3c 
can be safely ignored. 

Next we have fitted the excess energy per area at the bulk critical point with the 
ansatz 

Ee.(Lo, (3c) = B + a (Lo + 1,)-^+'/'' , (38) 

where we have used EhuikWc) = 0.602111(1) (Ref. |21|]) to compute Eex{LQ, Pc) and 
we have fixed z/ = 0.63002 (Ref. \12^)- The parameters of the fit are B, a and 
Lg. Note that B corresponds to a correction of the analytic background caused 
by the boundaries that only depends on the local properties of the system at the 
boundaries and therefore takes the same value for ++ and H — boundary conditions. 
In order to estimate errors due to subleading corrections we have also fitted with 

Ee.(Lo, (3c) = B + a (Lo + 1^)-^+'/" x (1 + c (Lq + L,)"^) , (39) 

where we have included quadratic corrections. 

For ++ boundary conditions we get with the ansatz ( 138|) for Lo,mm = 8 the 
results B = 7.1893(3), a = -8.045(1), L, = 1.915(2) and xVd-o.f. = 0.79. Using 
the ansatz ([39]) and Lo,™„ = 6 we get S = 7.1888(2), a = -8.042(1), L, = 1.905(1) 
and xVd.o.f. = 0.96. 

Instead, for H — boundary conditions we get using the ansatz ( 138|) for Lo,mm = 13 
the results B = 7.1947(4), a = -12.207(2), L, = 1.966(3) and xVd-o.f.' = 0.60. 
Using ansatz ( 15^ we get for LQ„iin = 8 the results B = 7.1864(5), a = —12.156(3), 
Lg = 1.830(6) and x^/d.o.f. = 0.53. The results of the two ansatze (I38f39p differ by 
several standard deviations, indicating that the systematical error due to corrections 
to scaling is clearly larger than the statistical one. Here we try to estimate this error 
from the difference between the results of the two ansatze f p8|39|) . Furthermore we 
have redone the fits above using shifted values for the input parameters EtuikiPc) and 
u to estimate the effect of their uncertainty on our results. In particular we find that 
by using u = 0.63012 instead of z/ = 0.63002 the values of our fitparameters shift 
considerably. E.g. for ++ boundary conditions and Lo,mm = 8 using ansatz f l35]) 
we get B = 7.1912(3), a = -8.040(1), L, = 1.909(2) and xVd-o.f. = 0.78. Taking 
into account the results of both ++ and H — boundary conditions we arrive at 

B = 7.189(6) (40) 

Ls = 1.9(1) (41) 

a++ = -8.04(1) (42) 

a+_ = -12.18(3) , (43) 

where we have taken the error mainly from the difference between the two different 
ansatze for the H — boundary conditions. We notice that the result obtained for 
Lg is fully consistent with that obtained from the analysis of the second moment 
correlation length above. 



B. The magnetization profile at the critical point 

In order to determine the constant Ls we have studied the magnetization at 
z = xq — {Lq + l)/2 = 0, i.e. in the middle of the film, for ++ boundary conditions. 
In the case of odd Lq we did use directly the value of the magnetization aX z = 0. 
In the case of even Lq we extrapolated the values of m at 2; = 3/2 and z = 1/2 to 
z = 0, assuming a quadratic dependence on z. For example for Lq = 24, 25, 26, 
28, 30, and 32 we get m\^^^ = 0.248488(6), 0.243670(4), 0.239111(4), 0.230695(4), 
0.223091(4), and 0.216181(4), respectively. 

Following eq. (^^, we have fitted our data with the ansatz 

m|^^o = C^(Lo + L,)-^/- (44) 

where Cm and Ls are the parameters of the fit. Note that /3/z/ = (1 + ri)/2 follows 
from scaling relations among the critical exponents. In our fits, we have fixed 
r] = 0.03627 (Ref. [l2|). In order to check for the effect of possible corrections, we 
have used in addition 

^L=o = Cm {Lo + LsY^/'' X (1 + c (Lo + L,)-2) . (45) 



Fitting with the ansatz f l44j) we find that the result for Ls is slowly decreasing with an 
increasing minimal thickness Lo^min that is included into the fit. For Lo,mm = 20 we 
find that x^/d.o.f. is still larger than two. For Lo,mm = 24 we get Cm = 1.34250(10), 
Ls = 1.937(4) and xVd-o.f.= 0.34. We have redone the fit with r] = 0.03637 instead 
of the central value rj = 0.03627. We find that the effect on Cm and Lg is much 
less than the statistical errors quoted above. Fitting with the ansatz fll^ we find 
for Lo,min = 16 the results Cm = 1.34171(17), L, = 1.867(12) and xVd.o.f.= 0.55. 
Also here we find that the error due to the uncertainty of rj is small compared with 
the statistical error quoted. Our results for Lg are in very good agreement with 
those obtained above. 

Finally in figure [1] we plot L^'^jjm{z) as a function of z/L^^eff using Ls = 1.9 
and rj = 0.03627 for ++ and H — boundary conditions. To this end we have used all 
thicknesses available with Lq > 16. The statistical errors are much smaller than the 
symbols that are used. For z/LQ^^ff ~ 0.4 the points fall nicely on unique curves 
for ++ and H — boundary conditions, respectively. For larger values of 2; a small 
scattering of the data can be observed. As the boundary is approached, this means 
z — )■ 1/2, the curves for ++ and H — boundary conditions fall on top of each other. 

C. Casimir force at the critical point 

We have computed 

A/(Lo, /?,) = /(Lo + 1/2, /?,) - /(Lo - 1/2, /?,) (46) 

using the algorithm discussed in ref. J39|. We have simulated ++ and H — boundary 
conditions on lattices of the thicknesses Lq = 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 
13.5, 15.5, 19.5, 23.5, 27.5, 31.5 and 39.5. For all these simulations, we have used 
L ^ 8Lq. We have checked that this is sufficient to avoid finite L corrections. These 




0,eff 



./3/v 



FIG. 1. We plot LQ^fr'm{z) as a function of z/Lo^^ff, where z = xq — {Lq + 1)/2 gives the 
distance from the middle of the film. The effective thickness of the film is -^o,e// = Lq + L^ 
using Ls = 1.9. For ++ and +— boundary conditions, data for films with Lq > 16 are 
used. 



simulations took in total about 10 month of CPU-time on one core of a Quad-Core 
AMD Opteron(tm) Processor 2378 running at 2.4 GHz. As update we have used 
the local heat bath algorithm. For lack of time and the still moderate amount of 
CPU time that was spent here, we made no effort to implement cluster updates or 
to implement the method using the multispin coding technique. 
We have fitted our data with the ansatze 

A/(Lo, Pc) = fbuikWc) - 0{O) (Lo + L,)-3 (47) 

and in order to check for the effect of subleading corrections to scaling 



A/(Lo, /3,) = f,^ik{/3,) - 6(0) (Lo + L,)-' x (1 + c (Lq + L 



(48) 



Fitting with the ansatz fHTj) we get for the ++ boundary conditions and Lo^min = 
11.5 the resuhs UuikiPc) = -0.0757368(3), ^(0) = -0.815(10), L, = 1.86(5) 
and x^/d.o.f.= 0.36. Using the ansatz (148|) and -Lo,mm = 6.5 we get the results 
fbuikWc) = -0.0757370(2), ^(0) = -0.824(5), L, = 1.91(4) and xVd.o.i.= 0.51. 

Fitting with the ansatz ( H7|) we get for the H — boundary conditions and LQ^rnin = 
11.5 the results fbrnkil^c) = -0.0757368(2), ^(0) = 5.617(16), L, = 1.930(13) 
and x^/d.o.f.= 1.11. Using the ansatz (H5]l and Lo^rnin = 6.5 we get the results 
fbuikWc) = -0.0757368(2), ^(0) = 5.610(14), L, = 1.912(17) and xVd.o.f.= 0.81. 



We notice that the results for ftuikWc) obtained from the two different boundary 
conditions are consistent. We conclude 

fbuiki(3c) = -0.0757368(4) . (49) 

Also the values for Ls obtained here are fully consistent with the estimate Lg = 
1.9(1) found above. As our result for the finite size scaling functions at the critical 
point of the bulk system we quote 

^++(0) = -0.820(15) (50) 

^^_(0) = 5.613(20) . (51) 

Also here we have checked that the uncertainty of Pc can be safely ignored. For a 
comparison of these results with previous ones given in the literature, see table IIIII 
below. 



VII. NUMERICAL RESULTS FOR THE CASIMIR FORCE IN A LARGE 
RANGE OF TEMPERATURES 



Here we compute the Casimir force using the method discussed by Hucht |40 . 
The details of the implementation are similar to [41|, where we have studied the 
thermodynamic Casimir force for films with free boundary conditions in the three 
dimensional XY universality class. 

We have simulated the model for both types of boundary conditions and the 
thicknesses Lq = 8, 9, 16, 17, 32 and 33 for a large number of /3- values in the neigh- 
borhood of the critical point. In tables ID and HTl we give the /3-values at which we 
have simulated and the statistics of our runs for the H — and the ++ boundary 
conditions, respectively. In the case of H — boundary conditions we also give the 
lattice size L that was used. Since for H — boundary conditions the correlation 
length is increasing with increasing /3 also L has to increase with increasing /3. In 
contrast, for ++ boundary conditions, the correlation length stays rather small for 
all temperatures. It has a maximum quite close to the critical point. Therefore we 
have used L = 32 for Lq = 8, 9, L = 64 for Lq = 16, 17 and L = 128 for Lq = 32, 
33 at all values of /3, where we have simulated at. 

We have measured the energy per area. Using these data we have computed 

AE{Lo, 13) = E{Lo + 1/2, /3) - E{Lo - 1/2, f3) - E^^iM . (52) 

The value for the energy density of the bulk system EbuikiP) is taken from sim- 
ulations of L^ or 2L X L^ lattices with periodic boundary conditions in all three 
directions. The linear lattice size L is taken sufficiently large to avoid significant 
finite size effects. For most values of /3 simulated here we have also a direct mea- 
surement of Ebuik{(3)- In a small neighborhood of /3c we have used instead the result 
of a fit with the ansatz 

E^ij,{P)=Ens + Cns{P-Pc) + a±\/3-/3c\'--'' + dns{/3-/3c)^ + b±\/3-Pc\''-'' . (53) 



For a discussion see section IV A of J2l|. Throughout the statistical error of Ehuik{(3) 
is clearly smaller than that oi E{Lq + 1/2, (3) — E{Lq — 1/2, (3). Also the systematical 
error caused by the interpolation with the ansatz (155|1 can be safely ignored here. 



In order to obtain Af^x we have numerically integrated AE^x using the trape- 
zoidal rule: 

n— 1 ^ 

- A/,,(/3„) ^ Yl 2 (^^+1 ~ f^'^ {AE,x{f3i+i) + AE,xm) (54) 

i=0 

where /3j are the values of /3 we have simulated at. They are ordered such that 
/3j+i > Pi for all i. The starting point of the integration (3q is chosen such that 
AEf^x{Po) = within the statistical error. 

The estimate obtained from the integration is affected by statistical and sys- 
tematical errors. The statistical one can be easily computed, since the AEex are 
obtained from independent simulations: 

e^-AfUPn)) = ^-^^-^^e'lAEUM] + ^-^^^^f^e'[AEUf3n)] 



n-l 



4 



2 



e'[AE,x{m (55) 



where e^ denotes the square of the statistical error. 

In order to estimate the error due to the finite step size /3j+i — /3j we have redone 
the integration, skipping every second value of (3; i.e. doubling the step size. We 
find that the finite step size errors are at most of the size of the statistical ones. 

In figures [2] and [3] we have plotted our results for the finite size scaling functions 
^+_(x) and ^++(x), respectively. The solid lines that are plotted linearly interpolate 
between the data points that we have computed. Note that the statistical error of 
Afex{Lo)LQ is of similar size as the thickness of the line. In both cases, in the upper 
figure we do not take into account any correction to scaling. This means we plot 
-Afex{Lo)Ll as a function of t[Lo/^o]^/^ using u = 0.63002. 

Not taking into account any correction, we see for both ++ and H — boundary 
conditions a clear discrepancy between the curves for Lq = 8.5, 16.5 and 32.5. 

Therefore in the lower part of the figures [2] and [3] we have replaced Lq by -^o,e// = 
Lq + Ls, using the value Lg = 1.9 obtained above from the finite size scaling study 
at the bulk critical point. This means that we have plotted —Afex{Lo){Lo + LsY as 
a function of t[(Lo + Ls)/^o]^'''^- Now the curves essentially fall on top of each other. 
Therefore we do not consider further corrections and take the curves obtained for 
Lq = 16.5 and 32.5 as our final result. The remaining small difference between 
Lq = 16.5 and 32.5 gives us some measure for the systematical error of our final 
result. 

Now let us discuss the properties of 6'++(a;) and 9+_{x). We see that 6'++ (a;) is 

negative and 6^ (x) is positive in the whole range of x. This means that in the 

case of ++ boundary conditions the force is attractive, while for H — boundary 
conditions it is repulsive. In both cases the function shows a single extremum. In 
the case of -|--|- boundary conditions it is located in the high temperature phase, 
while for H — it is in the low temperature phase. In order to accurately locate these 
extrema, we have computed the zeros of AE{Lq, P). For -|--|- boundary conditions 
we find /3^i„ = 0.37407(3), 0.38219(2) and 0.38569(2) for Lq = 8.5, 16.5 and 32.5, 
respectively. For these values of Pmin we have computed Xmin = ^mm[(-^o + -^s)/^o]^'''^ 
and correspondingly 9min = —Afex{(3min){LQ + Lg)^. As our final result we take the 



TABLE I. Statistics of our runs for the H — boundary conditions. In the first column we 
give the thickness that is considered, where for example Lq = 8.5 means that we have 
simulated films of the thicknesses Lq = 8 and 9. In the second column we give the linear 
extension L of the lattice in 1 and 2 direction. We have simulated at /3j = /3m,m + ^A/3 
in the interval [f3mim f^max]- In the last column we give the number of measurements for 
each of the simulations. 



Lq L jSmin 


Pmax 


A^ 


stat 


8.5 32 0.25 


0.325 


0.005 


200 000 


8.5 32 0.33 


0.348 


0.002 


200 000 


8.5 32 0.35 


0.38 


0.001 


200 000 


8.5 32 0.381 


0.385 


0.001 


300 000 


8.5 64 0.385 


0.43 


0.001 


150 000 


8.5 96 0.43 


0.46 


0.002 


100 000 


8.5 128 0.46 


0.5 


0.002 


100 000 


8.5 256 0.505 


0.56 


0.005 


100 000 


16.5 64 0.34 


0.348 


0.002 


200 000 


16.5 64 0.35 


0.384 


0.001 


200 000 


16.5 64 0.385 


0.395 


0.0005 


200 000 


16.5 128 0.395 


0.41 


0.001 


100 000 


16.5 256 0.412 


0.42 


0.002 


100 000 


16.5 512 0.422 


0.43 


0.002 


100 000 


16.5 512 0.44 


0.44 


0.01 


100 000 



32.5 128 0.36 0.355 0.005 1 000 000 

32.5 128 0.365 0.368 0.001 1 000 000 

32.5 128 0.369 0.3875 0.0005 1 000 000 

32.5 128 0.3875 0.39125 0.00025 1 000 000 

32.5 256 0.3915 0.395 0.0005 250 000 



value obtained for Lq = 32.5 using L^ 


= 1.9, u = 


0.63002 and ^o = 


= 0.2282. We 


arrive at 








x++,^i„ = 5.82(10) 


"++,min 


-1.76(3) , 


(56) 



where the quoted error takes into account the statistical error and the errors due 
to the uncertainties of Lg, ^o and u. 

For H — boundary conditions we find (3max = 0.39961(2), 0.39256(2) and 
0.389525(10) for Lq = 8.5, 16.5 and 32.5, respectively. In the same way as above 
for ++ boundary conditions we arrive at 

x+.,max = -5.17(7) 0+_,„,,. = 6.56(10) . (57) 

At the bulk critical point we get ^++(0) = 0.84(2) and ^+-(0) = 5.56(7). These 
results are less precise but fully consistent with those obtained in the previous 
section, eqs. O50f5ip . 
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FIG. 2. H — boundary conditions. In the upper part of the figure we plot — LgA/ea; as 
a function of t(Lo/^o)^ for Lq = 8.5, 16.5 and 32.5, where we use v = 0.63002 and 
^0 = 0.2282. In the lower part we have replaced Lq by L^^eff = Lq + Lg with Lg = 1.9. 
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FIG. 3. -|--|- boundary conditions. In the upper part of the figure we plot — LqA/ei: as 
a function of t(Lo/^o)^ for Lq = 8.5, 16.5 and 32.5, where we use v = 0.63002 and 
^0 = 0.2282. In the lower part we have replaced Lq by -^o,e// = Lq + L^ with Lg = 1.9. 
For a discussion see the text. 



TABLE II. Statistics of our runs for the ++ boundary conditions. The notation is the 
same as in the previous table for -\ — boundary conditions. Here we have used L = 
4(Lo - 1/2) for all values of 13. 



^0 Pmin 


Pmax 


A/3 


stat 


8.5 0.25 


0.295 


0.005 


5 000 000 


8.5 0.3 


0.348 


0.002 


5 000 000 


8.5 0.35 


0.358 


0.002 


10 000 000 


8.5 0.36 


0.378 


0.001 


10 000 000 


8.5 0.379 


0.395 


0.0005 


10 000 000 


8.5 0.396 


0.409 


0.001 


10 000 000 


8.5 0.41 


0.43 


0.002 


10 000 000 


16.5 0.31 


0.33 


0.01 


10 000 000 


16.5 0.34 


0.352 


0.002 


10 000 000 


16.5 0.354 


0.379 


0.001 


10 000 000 


16.5 0.38 


0.382 


0.0005 


10 000 000 



16.5 0.3825 0.39225 0.00025 10 000 000 
16.5 0.393 0.399 0.001 10 000 000 
16.5 0.4 0.406 0.002 10 000 000 



32.5 0.37 


0.375 


0.001 


10 000 000 


32.5 0.376 


0.3795 


0.0005 


10 000 000 


32.5 0.38 


0.3856 


0.0002 


10 000 000 


32.5 0.3858 0.3889 


0.0001 


10 000 000 


32.5 0.389 


0.3918 


0.0002 


10 000 000 


32.5 0.392 


0.3945 


0.0005 


10 000 000 


32.5 0.395 


0.396 


0.001 


10 000 000 



In ref. [42[ we have demonstrated at the example of films with periodic and 
free boundary conditions in the three dimensional XY universality class that the 
relation 9{x) = 2h{x) — ^h'{x), eq. (fT9l) . can be employed to compute 9{x) from the 
excess energy per area of the film, without taking the derivative with respect to the 
thickness Lq of the film. 

The main practical problem of this approach is that for free boundary conditions 
as well as symmetry breaking boundary conditions that are studied here, the ana- 
lytic part of the free energy per area and hence also of the energy per area suffers 
from a boundary correction that is not described by -^o,e// = Lq + Ls of the singular 
part. In section IVI Al we have already determined the value of this correction at 
the bulk critical point. However it turns out that it is not sufficient here to approx- 
imate this correction by a constant. Even by adding a term linear in the reduced 
temperature t to the analytic boundary correction, we could not reliably compute 
^++(x) and 6^_{x). We made no attempt to improve this by adding higher order 
terms. 



o 




FIG. 4. We plot our numerical results for ^-|_-|-(x) and —0-\ (x) obtained with Lq = 16.5 

and 32.5 using L^ = 1.9 for x > 0. For comparison we give the result of eq. ([3T]) . setting 
C^ = 1.5. For a discussion see the text. 



A. Behavior at large 



In figure m we have plotted 6^^{x) and —6^^{x) in the high temperature phase. 
For comparison we have plotted 6^^{x) given by eq. fl3T|) . We have fixed the constant 
C^ by matching the value at x ~ 20, where 0++(x) and —6^^{x) still agree within 
the error bars. We find 



C^ = 1.5(1) 



(5J 



Indeed for x ^ 20 at the level of our accuracy 6'++(x) and — ^+_(x) are equal. 
In the same range, the two curves are well approximated by eq. (1311) . 

Next let us turn to the low temperature phase. We have matched eq. ( 133|) with 
our numerical results obtained for Lq = 16.5 and 32.5 and ++ boundary conditions 
at a; fa —7. We get 



C^ = 0.20(5) . 



(59) 



As one can see from figure [5] there is reasonable match between our numerical results 
for 6'++(x) and eq. (!33l) for x ^ —5. In figure [5] we have plotted the statistical error 
of our results. The fact that for small x, within less than two standard deviations, 
the estimate of 9++{x) computed for Lq = 16.5 and Lq = 32.5 becomes equal to 
zero is a non-trivial validation of our numerical integration. 



CT> -0.2 - 




FIG. 5. We plot our numerical results for 0++(x) obtained with Lq = 16.5 and 32.5 using 
Ls = 1.9 for x < 0. For comparison we give the result of eq. (f33l) . setting C^ = 0.2. For 
a discussion see the text. 



B. Correlation length of the films 



For all simulations discussed above we have measured the second moment cor- 
relation length as defined in section IIIICI The correlation length is interesting for 
practical purpose, since we have to choose the lattice size L in 1 and 2 direction 
such that L ^ E,2nd in order to avoid sizable effectively two dimensional finite size 
effects. Furthermore we shall discuss the finite size scaling behavior of the second 
moment correlation length of the film to further probe the theoretical expectations 
on corrections to scaling. 

To this end, we have plotted in figure [6] for ++ boundary conditions ^2nd/-^o,e// 
of the film as a function of the scaling variable x = t[Lo,e///Co]^'''' for the thicknesses 
Lq = 8, 9, 16, 17, 32 and 33. Using Lg = 1.9 instead of L^ = clearly improves the 
collapse of the curves obtained from different thicknesses Lq. Using Lg = 1.9, in 
the range —20 ^ x ^ 20 the curves obtained for different thicknesses fall on top of 
each other within the error bars. For larger values of x there is some discrepancy 
between the thicknesses Lq = 8 and 9 and Lq = 16, 17, 32 and 33 on the other 
hand. This can be attributed to analytic corrections to scaling. For all thicknesses 
C,2nd/LQ^eff assumes a single maximum at x ~ 7. 

Figure [7] is the analogue of figure [6] for H — instead of ++ boundary conditions. 
Also here we find, using Lg = 1.9 a nice collapse of the curves obtained for the 
different thicknesses of the films. Now ^2nd/-^o,e// is monotonically increasing with 
decreasing x. In figure [7] we have stopped, a bit arbitrary, at a; = —50. For 
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FIG. 6. For ++ boundary conditions, we plot C2nd/-^o,e// as a function of the scaling 
variable x = t[Lo,e///'^o]^ for Lq = 8, 9, 16, 17, 32 and 33 using Lg = 1.9. Notice 
that ^2nd is the second moment correlation length of the film, while ^o appearing in the 
scaling variable x is the amplitude of the correlation length of the bulk system in the high 
temperature phase. 



X ~ —79.7, the smallest value of x that we have reached for Lq = 9, we get 

i2nd/ Lo^eff ~ 3.5. 

With an increasing correlation length the autocorrelation time of the Metropo- 
lis update increases. Therefore simulations become increasingly difficult as we go 
deeper into the low temperature phase, towards smaller values of x. As a conse- 
quence we had to stop at x ~ —20.4 and —21.3 for Lq = 32 and 33, respectively. 

VIII. COMPARISON WITH OTHER THEORETICAL RESULTS AND 
EXPERIMENTS 



The scaling functions 6'++ and ^+_ have been computed recently by using Monte 



Carlo simulations of the spin- 1/2 Ising model on the simple cubic lattice 16, 17 
The results are presented in figures 3 and 4 of [16| and 9 and 10 of [17[ for -|--|- and 
H — boundary conditions, respectively. For both types of boundary conditions, the 
final result depends strongly on the precise form of the ansatz, see eqs. (18,20,21,23) 
of [iTJ, for corrections to scaling that is chosen. Qualitatively, the curves for both 
-|— |- and H — boundary conditions agree with ours. For the position of the extrema 
the authors of [iTJ quote x++^mm = 5.90(8) and x+_^max = ~5.4(1) in the caption 
of their figures 9 and 10, respectively. These are in quite good agreement with our 
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FIG. 7. Same as figure [6] for H — instead of -|--|- boundary conditions. 



results. In J4J, |45|, see the discussion below eq. (14) of [4q], the authors extract 



the amplitude C^ from the data of [17|]. Their result depends on the ansatz that is 
chosen for the corrections and also on the boundary conditions. Using the ansatz 
that is denoted by (i) in figures 9 and 10 of Q, they find C^ = 1.51(2) and 1.82(2) 
for ++ and H — boundary conditions, respectively. Instead, using the ansatz that 
is denoted by (ii) they arrive at C^ = 1.16(2) and 1.38(2), respectively. It is clear 
from these numbers that systematical errors due to corrections to scaling are much 
larger than statistical errors. Taking this into account, there is nice agreement with 
our estimate C^ = 1.5(1), eq. fl58|) . 

In figure [8] we compare our result for 6'++(x) with that obtained by using the 
de Gennes-Fisher local- functional method [19(. As input the method uses universal 
amplitude ratios of the bulk system. Here we made no effort to redo the calculations 
of [l9| using our updated values for the universal amplitude ratios [21j] and value for 
the exponent u (Ref. [l2|). Instead, we have copied the curve from figure 1 of 19 . 
Overall we find a reasonable agreement with our result. We see a very small shift of 
the local-functional method curve towards larger values of x compared with ours. 
Clearly, the value of the minimum of the curve obtained by the local-functional 
method is smaller than that of ours. 

The authors of |43| have studied wetting films of a binary mixture of methylcyclo- 
hexane and perfluoromethylcyclohexane. They have deduced the thermodynamic 
Casimir force from measurements of the thickness of the film. Their result for 
^+_ (x) given in figure 3 of |43| is more or less consistent with but much less precise 
than our result. The authors of 4J, |45| have studied the thermodynamic Casimir 
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FIG. 8. We plot the result of ref. (I9j] for ^_|__|_(x) obtained by using the de Gennes- 
Fisher local-functional (LF) method. We have copied the curve from fig. 1 of 19|]. For 
comparison we plot our numerical results for ^^^(x) obtained with Lq = 16.5 and 32.5 
using Ls = 1.9. 



force between colloidal particles that are immersed into a mixture of water and 
lutidine and the surface of the cell. The surface of the particle was prepared such 
that it either preferentially absorbs water or lutidine. Hence both ++ and H — 
boundary conditions were accessible. A major problem in the interpretation of the 
experimental data is to disentangle the thermodynamic Casimir force from other 
forces. It turns out that only for relatively large x, reliable results could be ob- 
tained. Theoretically the colloid al p article and the surface of the cell are described 
by a sphere and a plane. In J4^ |45| the Derjaguin approximation had been used to 



obtain a prediction for this geometry starting from the theoretical results for the 
universal finite size scaling functions 6^^{x) and 6^_{x) for the film geometry. The 
45| have fitted their data with the equivalent of ansatz (EI]), taking 



authors of |4 

^0 as free parameter. Their result for ^o is consistent with that obtained from the 
analysis of bulk quantities. This check could be made more stringent by replacing 
the theoretical estimate of C^ of [4^, |45| by ours eq. ( 158|) . 

Finally in table IIIII we have summarized results obtained for the scaling func- 
tions at the bulk critical point. In the literature, results obtained by field theoretic 
methods l32|, the de Gennes-Fisher local- functional method [l8|, Monte Carlo sim- 
ulations I id . llTI . |32| and experiment |43| can be found. Mostly, in the original 
work, the so called Casimir amplitude A = 9{0)/2 is quoted. We see that field 
theoretic methods, in particular the e-expansion, are not able to provide quanti- 
tatively satisfying results. Those of the de Gennes-Fisher local-functional method 



TABLE III. Comparison of our results for 9-^-^(0) and 9j^ (0) with those given in the 

hterature. For a discussion see the text. 



Ref. Method ^++(0) ^+-(0) 



[321 e-expansion -0.346 3.16 
[32] d = 3 expansion -0.652 4.78 
[18] local-functional -0.84(16) 6.2 



[43J 


experiment 


6(2) 


[32j 


Monte Carlo 


-0.690(32) 4.900(64) 


[16J 


Monte Carlo 


-0.884(16) 5.97(2) 


m 


Monte Carlo 


-0.75(6) 5.42(4) 


here 


Monte Carlo 


-0.820(15) 5.613(20) 



181 ] are in much better agreement with ours. The results of previous Monte Carlo 



simulations differ by more than the quoted error bars from our results. Note that 



are 



in ref. ^6[ only the statistical error is quoted. The numbers quoted for p/7 
obtained by using an ansatz different from that of [16|, which explains the difference 
between them. In figure 8 of [l7| the authors give in addition to the results obtained 
with their preferred ansatz those obtained by using two alternative ansatze. From 
this comparison one might conclude that the systematical error is larger than the 
statistical one that we quote in table lllli 

As we have seen here, for the thicknesses that can be studied today, corrections 
to scaling, in particular those caused by the boundaries, are numerically important. 
In order to get an accurate result for the scaling limit, these corrections have to 
be properly taken into account. In the generic case, when corrections oc Lq'^ , with 
u = 0.832(6), and oc L^^ are present this is a difficult task. 



IX. SUMMARY AND CONCLUSIONS 

We have studied the thermodynamic Casimir force for thin films in the three di- 
mensional Ising universality class. In particular we have studied symmetry breaking 
boundary conditions. We consider the two cases ++ and H — , where the fixed spins 
at the boundary are either all positive or are positive at one boundary and negative 
at the other. We have simulated the improved Blume-Capel model on the simple 
cubic lattice. The boundary conditions are expected to cause corrections that are 
to leading order oc Lq^. In general it is hard to disentangle such corrections from 
leading corrections to finite size scaling which are oc Lq'^ where u = 0.832(6) (Ref. 
[12|). In the improved model, corrections to scaling oc Lq^ are eliminated. This fact 
very much simplifies the analysis of the Monte Carlo data. In particular we could 
clearly demonstrate that the corrections caused by the boundaries can be expressed 
by an effective thickness -^o,e// = -^o + Lg. For our model we find, for both ++ and 
H — boundary conditions Lg = 1.9(1). 

Having corrections to scaling well under control, we have obtained reliable re- 
sults for the universal finite size scaling functions 6'++(x) and 6'+_(a;), where x = 



t[Lo,e///'Co]^^'^) of the thermodynamic Casimir force. For large values of x, we have 
compared our estimates for ^++(x) and 9+-{x) with the prediction fl3T|) derived by 
using the transfer matrix formalism. We find good agreement. For large values 
of —X we have compared 9++{x) with eq. ([33D also derived by using the transfer 
matrix formalism. Also here we find agreement. 

Finally we have compared our estimates for 9++{x) and 9^^{x) with field the- 
oretic calculations, the de Gennes-Fisher local-field method, previous Monte Carlo 
simulations and experiments. While field theory does not provide quantitatively 
satisfying results, those of the local-field method are in quite reasonable agreement 
with ours. Also the results of previous Monte Carlo simulations are essentially in 
agreement with ours. 
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